New method for the 3D Ising model 
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A simple, general and practically exact method is developed for the equilibrium properties of the 
macroscopic physical systems with translational symmetry. Applied to the Ising model in two and 
three dimension, a modest calculation gives the spontaneous magnetization and the specific heat to 
less than 1% error. 
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To calculate the equilibrium properties of a given phys- 
ical system will probably be the most basic theoreti- 
cal task. In 1944, Onsager [l| calculated the partition 
function of the two dimensional Ising model analytically 
and demonstrated the power of the exact solution, fol- 
lowed by the exact spontaneous magnetization by Yang 
0. Since then some remarkable progresses have been 
made in two-dimensional classical and one-dimensional 
quantum systems analytically and by numerical 

methods such as DMRG (density matrix renormalization 
group) [7|. However, reflecting the three-dimensionality 
of our world, much needed knowledge is for three dimen- 
sion, and several thousand papers have been said written 
on the three dimensional Ising model. Among others, a 
tremendous effort has been paid to the precise determi- 
nation of the critical temperature T c in association with 
the concept of universality near the critical point [H . Ow- 
ing to the concerted effort of RG-finite size scaling- Monte 
Carlo (MC) method, T c is known as precise as 4.5114(1) 
@, as well as the spontaneous magnetization |lo| . The 
effort of high and low temperature expansions by the di- 
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agrammatic method for T c is also remarkable 
another noteworthy recent progress is along the line of 
DMRG El. 



We here present a new method for the equilibrium 
properties of the macroscopic physical systems with 
translational symmetry. Our method is purely algebraic, 
works directly on a oo x oo x oo lattice and, other than 
seeking a convergence in entanglement space (see below) , 
it does not invoke any other notions such as numerical 
RG, nor make any approximations. We have calculated 
the spontaneous magnetization and the specific heat in 
two and three dimensional Ising models to less than 1% 
error. Note that without the help of RG, the 1% preci- 
sion is what one can typically expect from exact methods 
such as Bethe Ansatz 0, The new method is simple, 
general and computationally efficient. In fact, our 30 
minutes calculation already gives fairly precise results as 
demonstrated below. 



square or simple cubic lattice, 
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where ai takes + or - 1 and the summation is over the 
nearest neighbor pairs. The partition function is given by 
Z = Tr[exp(— /3H)], where Tr means to take trace over 
the 2^ spin configurations with N being the number of 
lattice sites and /3 — 1/kT. We follow the following steps 
for the two dimensional case. 

First, note that the local pair-density matrix can be 
written as exp(/3criC7j) = cosh(/3) + sinh(/3)<7j<7j. This 
is the simplest case of more general statement that any 
local density matrix when regarded as a real symmet- 
ric matrix can be written by singular value decomposi- 
tion (SVD) as Aij = J2k Vik^kV jk , where V k and A fc are 
eigenvectors and eigenvalues of the local density matrix 
A. Denoting these two terms by a "bond index" i = 1,2, 
the partition function is written as a summation over all 
the bond index configurations on the square lattice. Now 
each term in this summation is a product over all the 
sites of a local spin operator constructed from Ai and 
\2CTi, depending on the four bond indexes surrounding 
the i-th spin, where Ai = ^/cosh(/3) and A2 = ^/sinh(/3). 
For example at site i, if the surrounding four bond in- 
dexes are all 1, then the local spin operator at the site 
is Af and its trace is, Tr\\ — 2\\. If one of the four 
bond indexes is replaced by 2, then the local spin trace 
is, Tr\\ \ 2 <Ji = 0. Let us denote these local spin traces 
as Tijki with the four bond indexes i, j, k and I taking 1 
or 2. Clearly, T ljM = 2Af , 2A?A^, 2A| for i + j + k + I = 
4,6, and 8 respectively and zero otherwise. Performing 
the spin trace over all the sites, we have 
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Second, by virtue of infinite system, we can conveniently 
assume the periodic boundary condition in both hori- 
zontal and vertical directions. Let us denote the leftmost 
horizontal bond indexes by a vector ai and the right- 
most ones by aN+i- The periodic boundary condition is 



Let us consider the ferromagnetic Ising model on a then expressed by the Croncckcr delta S : 



aia N + i 



which can 
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be expanded by an orthonormal complete set {<p n (a)} 

as 5 aiaN+1 = J^nVni^Vni^+l)- N ° te that a11 tlie 

quantities in this paper are real. Let us now denote the 
rightmost column composed of a vertical one-dimensional 
array of Ts as K (ajv, a\+i). Note that we can certainly 
regard this quantity as a function of the two bond-index 
vectors ajy and aN+i after summing over the vertical 
bond indexes associated with this column. Now choose 
the above introduced eigenstate as the one which satisfies 
a transfer-matrix eigenvalue equation 

iiT(a N ,a N+ i)^, i (a N+ i) = /x„<£„(a N ) (3) 

where and below the repeated indexes imply a sum- 
mation. Repeatedly using Eq(|3]), the partition func- 
tion is evaluated as Z = J2 n f-n ~ * Po i R tne large N 
limit, where the suffix indicates the largest eigenvalue. 
The specific heat can be obtained from the free energy 
F = —kTln(Z) by differentiating twice with respect to 
temperature. To avoid unnecessary error due to numer- 
ical differentiation, however, one might calculate the in- 
ternal energy per site, E = —D < Oi<jj > where D is the 
space dimensionality and < OiOj > means the correla- 
tion of nearest neighbor spins, and then the specific heat 
is obtained by differentiating the energy with respect to 
temperature once. Likewise the spontaneous magnetiza- 
tion per site is given by M =< a, >. Third, we write 
<^o( a N+i) in the same form as the operator K(a^, a^+i) 
with each r^jy being replaced by a three index quantity 
Cimm thus without the index aj\j- This is justified as a 
Schmidt decomposition [3] of eigen functions with trans- 
lational invariance. Or more simply, one can derive this 
form by a successive use of SVD. Consider, for example, 
a wave function (^(01020304). Regarding this as a ma- 
trix of the left index a± and the right index {020304}, 
we have ^(01020304) = J2 a A a ia p a B {a2aaai}a . The 
quantity B can in turn be regarded as a matrix of 
the left index {02a} and the right index {0304}, thus 
SVD gives B{oa0304}e« = J2 l 3 C {a 2 a}p^l3 D {a 3 a i }i3- Like- 
wise, D {a3ai}l3 = E 7 £, {a 3 /3}7 A 7 i? a 4 7- Putting together, 
rewriting A aia as A a (a 1 ), C {a2Q yfj as C a( a(a 2 ), £{a 3/ 3} 7 
as Efj 1 {a3) and F ail as F 7 (a4), and appropriately ab- 
sorbing p a , \[j and A 7 into the matrixes A, C, E and F, 
one gets 93(01020304) = A a (ai)C a i3(a 2 )E l 3- f (a 3 )F 7 (a 4: ). 
For an infinite system with translational symmetry and 
periodic boundary condition, one arrives at the claimed 
form, with a, j3 and Oj-i, e.g., corresponding to I, m, 
and n in the above Q mn , 

Vo(aN+i) = ■■.(af3(ai-i)(f3~ / {a i )(~ f 5(a i+1 )... (4) 

Fig.l represents the transfer matrix eigenvalue equation 
EqQ schematically. In quantum information theory, 
these indexes a, j3 and 7 are known as entanglement 
[lij . Considering only 1 for these indexes is a simple 
mean- field-like approximation for v?o( a N+i)- Allowing 
larger values, one takes into account the effect of correla- 
tion with increasing precision. An important note here is 
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FIG. 1: Schematic figure of the transfer matrix eigenvalue 
equation, EqQ. 

that Eq(|3]) is not peculiar to the Ising models, but rather 
a general statement for any macroscopic systems with 
translational symmetry. Another important note is, as 
is evident from the derivation above, the range of values 
which the entanglements a, f3 and 7 can take, grows ex- 
ponentially for increasingly large systems. Thus the suc- 
cess of the present method crucially relies on a rapid con- 
vergence of physical quantities with the increase of entan- 
glement. Fourth, we handle the eigenvalue problem (3) as 
a variational problem, namely we maximize the quantity 
po = (poKipo/ifofo by iteration starting with an input 
state for tpo. First consider the numerator. A local ingre- 
dient of this quantity is A U > e , mm 'f = Clmn^nefn'Q'm'n'- 

Since A is real symmetric, we can use SVD again and 
write as A — XvX tr where the matrixes X and v are 
made up of eigenvectors and eigenvalues of A and tr 
means the transpose. Since X is orthogonal, the summa- 
tion over the combined entanglement-bond indexes We 
in the numerator can be done for N-f times, and thus we 
only keep the largest eigenvalue vq and eigenstate Xo- We 
have (poKtpo — VQ~ 1 xo tr Axo. The denominator can be 
handled likewise. Let us denote the corresponding largest 
eigenvalue and eigenstate as po and yo- Note that po 
now contanins Qmn in quadratic form. Maximizing this 
quantity with respect to Q mn then leads to a generalized 
eigenvalue problem: 

%0W e-T enn' /%mm' fCl'm'n' ~ PoVOW Snri'VOmm' 0' m' n' (5) 

where po = p,Qi> ~ x /p$~ . We solve Eq(5) for the next 
Q and continue until convergence. Finally after the con- 
vergence, we can calculate the internal energy and spon- 
taneous magnetization in terms of the obtained eigen- 
values and eigenstates. The procedure is the same as 
above with only difference being one (two) of Ts in the 
expression for the partition function should be replaced 
by Tijki = 2\\ X 2 , 2A1A2 for i+j+k+l = 5, 7 respectively 
and zero otherwise for calculating < <x; > (< o~i(Tj >). 
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We have, 

M = x tr ix /(AoPo) (6) 
where A is the same as A with T replaced by T. Likewise 

E = - Wi 2 x /(MoPo^o) (7) 

The extension to three dimension is straightfor- 
ward, r and r now have 6 indexes, Tijki mn = 
2Af,2A|A^2AfA^,2A| for i +j + A; + I + m + n = 6,8,10 
and 12 respectively and zero otherwise and Tijkimn = 
2Af A 2 , 2A?A^, 2AiA| for i + j + k + I + m + n = 7,9 and 
11 respectively and zero otherwise, and the bond index 
vector a represents a two dimensional array. The largest 
eigenstate ipo (a) of the operator K is now expressed as a 
product of Cabcda j the first four indexes denoting entan- 
glement while the last one a bond index. For notational 
simplicity, the suffix zero for denoting the largest eigen- 
value and eigenstate will be dropped from now on. The 
maximization of fi = (pKip/iptp now takes the structure 
of a " Russian dole" . In fact the calculation of the numer- 
ator is identical with that of the partition function in two 
dimension with T replaced by B — Cabcda^aefghpCijkip- 
The bond indexes are replaced by the composite of bond- 
entanglement indexes (aei), (bfj) etc. which we will de- 
note by a single index below when possible. B can then 
be expressed like -Bijjy. We consider an eigenvalue prob- 
lem similar to EqQ. Let us write the largest eigenvalue 
as A and the local constituent of the eigenstate as r) tt 'i 
(corresponds to Q m n in two dimension) where the first 
two indexes represent the new entanglement indexes and 
the last a composite of bond-entanglement indexes. Let 
us denote the quantities corresponding to uo, Xo and po 
in two dimension as r, u and e. The denominator is like- 
wise treated. Here L is replaced by C = CabcdaCijkia, 
and the composite of bond-entanglement indexes is now 
(ai) etc. Corresponding to A,tj,t, u and e, we write 
f2, £,Td,v and e^. It is convenient to define a quantity 
G p u,p't'l = Vpp'jBijki'qtt'k and associated G which is ob- 
tained from G with F replaced by F. Likewise define 
O p u,p'tn = ipp'jCijkiitt'k- With these notations, we have 

fi = Wu tr Gu/v tr Ov (8) 

with W = A N ~ 1 T N ~ 1 e^/(e N n N - 1 T^- 1 ). Note that C is 
quadratically involved here, and therefore the maximiza- 
tion of /i gives a similar generalized eigenvalue problem 
as Eq©, 

UptqsrTJpp' abc^- arr' cc' fiVtt' a'b' c'^p' t' q' s' r' Csbs'b' @ 

= iiVptqs^pp'abSa(l£,tt'a'b'Vp't'q's'Csbs'b'l3 (9) 

where fi = pW . When the iteration is converged, the 
spontaneous magnetization and the internal energy in 
three dimension are calculated as, note the same formula 
EqQ and Eq((7|) in two dimension, 

M = u tr Gu/(pr d ) (10) 



— r -*1— <,l r 


1 1 1 1 1 

\ 


- 


1 . 1 


, 1,1, 


\ 

■ 

J 
I 

• 1 ■ 1 1 ■ 1 



1.0 1.5 2.0 2.5 3.0 3.5 
Temperature 

FIG. 2: Spontaneous magnetization vs. temperature in two 
dimension. From right to left, the entanglement p=l, 2 and 8 
(star, red online). The leftmost solid line (blue online) is the 
exact result of Yang Q] . 




Temperature 

FIG. 3: The same as Fig.2 for the specific heat. The peak is 
at T = 2.278. The thick solid line (blue online) is the exact 
result of Onsager [l[ . 

E = -Du tr G 2 u/(pTT d ) (11) 

The results for the spontaneous magnetization and the 
specific heat in two dimension are shown in Fig.2 and 3. 
The results with the entanglement p = 8 are indistin- 
guishable from the exact ones of Onsager and Yang. The 
obtained critical temperature 2.278 is 0.4% off the exact 
value 2.269185. The results in three dimension are shown 
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FIG. 4: Spontaneous magnetization vs. temperature in three 
dimension. From right to left, the entanglement (p,q)=(l,l), 
(1,2), (2,1), (2,2) and (2,4) (star, red online), where p refers 
to the entanglement associated with £ and q with £ and n. 
For p=l, the (1,2) result is the converged one. The leftmost 
solid line (blue online) is the empirical formula from the RG- 
finite-size-scaling-MC method [Hj. 
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FIG. 5: The same as Fig. 4 for the specific heat. The peak is 
at T=4.547. 



in Fig. 4 and 5. The result for the spontaneous magneti- 
zation with the entanglement (p,q) = (2,4) is very close to 
the empirical formula due to the RG-finite-size-scaling- 
MC method 0, 

M(t) = -i°' 32694109 (1.6919045 
-0.34357731t - 50842026 - 0.42572366i), (12) 



for 0.0005 < t < 0.26 where t=l-0.2216544/fcT. The ob- 
tained critical temperature 4.547 is 0.8% off the believed 
exact value 4.5114(1). The specific heat in three dimen- 
sion has a marked difference from two dimension, namely 
it is highly asymmetric near the critical point. In fact 
our result with the entanglement (2,4) has quite a resem- 
blance to both the MC prediction [14| (also on page 479 
of [llj]) and the specific heat experiment on argon at the 
gas-liquid phase transition [3] (also on page 13 of 
It may be worth pointing out that the method can calcu- 
late equilibrium properties easily to a high precision. As 
is seen in Fig. 1-4, the p=2 case in two dimension and the 
(2,2) case in three dimension give fairly precise results, 
but it took only 30 minutes using a single PC of about 
1GHz processing speed. 

Because of simplicity and generality of the method, 
many applications and further progresses are anticipated, 
but from a theoretical viewpoint, the entanglement p = 
3 case in three dimension needs to be resolved first to see 
the swift convergence with the increase of entanglement. 
The simplest case (3,1), however, was done easily giving 
at most (near the transition point) 10~ 4 relative correc- 
tion to the case (2,1), indicating that the correction from 
p = 3 will indeed be very small. After confirming this 
point, the method would serve as an exact method for 
most of the purposes in equilibrium statistical physics. 
As for a much more accurate analysis of the transition 
point such as T c and various critical exponents beyond 
1% level, as was done in [§-[ll|, Qj}, it is interesting to 
see how a combination of RG and the present method 
performs. 

In conclusion, the essence of the new method shall be 
summarized and further clarified. First, it is not a kind of 
cluster mean field theories or transfer-matrix mean field 
theories as used, e.g., in the coherent anomaly method 
18]. In our word, the cluster mean field theories are de- 
scribed as the case with entanglement = 1, but each £ 
there represents a cluster of bonds. By considering entan- 
glement greater than 1, we go beyond mean field theories. 
Second, we have not made any assumptions, nor used any 
peculiar procedures only applicable to the Ising models. 
We have used SVD repeatedly and the key idea of our 
method is to fully implement the translational symmetry. 
Third, one might still press for an evidence of power of 
the new method. For that matter, we simply note that 
the new method has been successfully extended to the ID 
Hubbard model, reproducing the main results of Bethe 
Ansatz [l9[ • Finally, a deep question is if our method of- 
fers a new physical concept or picture. For this, we point 
out that the state ip® for the Hubbard model includes, as 
a special case, the resonating-valence-bond (RVB) state 
proposed by Anderson for the high-temperature super- 
conductors [2o( j . We note that the resonating character 
of the Anderson wave function is a natural consequence of 
the translational symmetry. In this sense, our state <po is 
always resonating. It could be any states, superconduct- 
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ing, insulating, matallic or magnetic. Our method thus 
opens a possibility of precise determination of the phase 
diagram of the 2D Hubbard model. An effort toward 
this goal is currently under way, and will be reported in 
a future publication. 
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